<!DOCTYPE html>
<html>
<head><meta charset="utf-8" />
<title>StokesP1bP1femrate</title><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.1.10/require.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/2.0.3/jquery.min.js"></script>

<style type="text/css">
/* Overrides of notebook CSS for static HTML export */
body {
  overflow: visible;
  padding: 8px;
}
div#notebook {
  overflow: visible;
  border-top: none;
}@media print {
  div.cell {
    display: block;
    page-break-inside: avoid;
  } 
  div.output_wrapper { 
    display: block;
    page-break-inside: avoid; 
  }
  div.output { 
    display: block;
    page-break-inside: avoid; 
  }
}
</style>

<!-- Custom stylesheet, it must be in the parent directory as the html file -->
<link rel="stylesheet" type="text/css" media="all" href="../doc.css" />
<link rel="stylesheet" type="text/css" media="all" href="doc.css" />

<!-- Loading mathjax macro -->
<!-- Load mathjax -->
    <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS_HTML"></script>
    <!-- MathJax configuration -->
    <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
        tex2jax: {
            inlineMath: [ ['$','$'], ["\\(","\\)"] ],
            displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
            processEscapes: true,
            processEnvironments: true
        },
        // Center justify equations in code and markdown cells. Elsewhere
        // we use CSS to left justify single line equations in code cells.
        displayAlign: 'center',
        "HTML-CSS": {
            styles: {'.MathJax_Display': {"margin": 0}},
            linebreaks: { automatic: true }
        }
    });
    </script>
    <!-- End of mathjax configuration --></head>
<body>
  <div tabindex="-1" id="notebook" class="border-box-sizing">
    <div class="container" id="notebook-container">

<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="P1bP1-Finite-Elements-for-Stokes-Equations">P1bP1 Finite Elements for Stokes Equations<a class="anchor-link" href="#P1bP1-Finite-Elements-for-Stokes-Equations">&#182;</a></h1>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>This example is to show the convergence of P1b-P1 finite elements for the Stokes equation on the unit square:</p>
$$- \Delta u + {\rm grad}\, p  = f \quad {\rm div}\, u    = 0  \quad  \text{ in } \quad \Omega,$$<p>with the pure Dirichlet boundary condition. The solver is based on a DGS type smoother.</p>
<p><strong>References</strong>:</p>
<ul>
<li><a href="http://www.math.uci.edu/~chenlong/226/FEMStokes.pdf">Finite Element Methods For Stokes Equations</a></li>
<li><a href="../project/projectFEM.html">Project: Finite Element Methods for Stokes Equations</a></li>
</ul>
<p><strong>Subroutines</strong>:</p>

<pre><code>- StokesP1bP1
- squareStokes
- femStokes
- Stokesfemrate

</code></pre>
<p>The method is implemented in <code>StokesP1bP1</code> subroutine and can be tested in <code>squareStokes</code>. Together with other elements (CRP0, P2P0, P2P1, isoP2P0, isoP2P1, P1bP1), <code>femStokes</code> provides a concise interface to solve Stokes equation. The P1b-P1 element is tested in <code>Stokesfemrate</code>.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="P1b-P1-element">P1b-P1 element<a class="anchor-link" href="#P1b-P1-element">&#182;</a></h2><p>The velocity is P1 Lagrange element enriched by a cubic bubble function and the pressure is P1 Lagrange element. This pair is also known as MINI element. The bubble function is introduced to stabilize the non-stable P1-P1 pair.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[1]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">clear</span> <span class="n">all</span><span class="p">;</span>
<span class="n">imatlab_export_fig</span><span class="p">(</span><span class="s">&#39;print-png&#39;</span><span class="p">)</span>  <span class="c">% Static png figures.</span>
<span class="c">%% Local indexing of DOFs</span>
<span class="n">node</span> <span class="p">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">;</span> <span class="mf">0.5</span><span class="p">,</span> <span class="nb">sqrt</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span><span class="o">/</span><span class="mi">2</span><span class="p">];</span>
<span class="n">elem</span> <span class="p">=</span> <span class="p">[</span><span class="mi">1</span> <span class="mi">2</span> <span class="mi">3</span><span class="p">];</span>
<span class="c">% elem2dof = 1:6;</span>
<span class="n">set</span><span class="p">(</span><span class="n">gcf</span><span class="p">,</span><span class="s">&#39;Units&#39;</span><span class="p">,</span><span class="s">&#39;normal&#39;</span><span class="p">);</span> 
<span class="n">set</span><span class="p">(</span><span class="n">gcf</span><span class="p">,</span><span class="s">&#39;Position&#39;</span><span class="p">,[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mf">0.3</span><span class="p">,</span><span class="mf">0.2</span><span class="p">]);</span>
<span class="n">subplot</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">showmesh</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span> <span class="n">findnode</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="s">&#39;all&#39;</span><span class="p">,</span><span class="s">&#39;noindex&#39;</span><span class="p">);</span> <span class="n">findelem</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">subplot</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">)</span>
<span class="n">showmesh</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span> <span class="n">findnode</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="s">&#39;all&#39;</span><span class="p">,</span><span class="s">&#39;noindex&#39;</span><span class="p">);</span>
</pre></div>

    </div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Dirichlet-boundary-condition">Dirichlet boundary condition<a class="anchor-link" href="#Dirichlet-boundary-condition">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[6]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="c">%% Setting</span>
<span class="c">% mesh</span>
<span class="p">[</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">]</span> <span class="p">=</span> <span class="n">squaremesh</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span><span class="mf">0.25</span><span class="p">);</span>
<span class="p">[</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">]</span> <span class="p">=</span> <span class="n">uniformrefine</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">);</span>
<span class="n">bdFlag</span> <span class="p">=</span> <span class="n">setboundary</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;Dirichlet&#39;</span><span class="p">);</span>
<span class="n">mesh</span> <span class="p">=</span> <span class="n">struct</span><span class="p">(</span><span class="s">&#39;node&#39;</span><span class="p">,</span><span class="n">node</span><span class="p">,</span><span class="s">&#39;elem&#39;</span><span class="p">,</span><span class="n">elem</span><span class="p">,</span><span class="s">&#39;bdFlag&#39;</span><span class="p">,</span><span class="n">bdFlag</span><span class="p">);</span>
<span class="c">% pde</span>
<span class="n">pde</span> <span class="p">=</span> <span class="n">Stokesdata1</span><span class="p">;</span> 
<span class="c">% options</span>
<span class="n">option</span><span class="p">.</span><span class="n">L0</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">maxIt</span> <span class="p">=</span> <span class="mi">4</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">printlevel</span> <span class="p">=</span> <span class="mi">1</span><span class="p">;</span>
<span class="n">option</span><span class="p">.</span><span class="n">solver</span> <span class="p">=</span> <span class="s">&#39;asmg&#39;</span><span class="p">;</span> <span class="c">% use auxiliary space preconditioner</span>
</pre></div>

    </div>
</div>
</div>

</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[7]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="n">option</span><span class="p">.</span><span class="n">elemType</span> <span class="p">=</span> <span class="s">&#39;P1bP1&#39;</span><span class="p">;</span>
<span class="n">femStokes</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">option</span><span class="p">);</span>
</pre></div>

    </div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

    <div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>#dof:   1763,  #nnz:  12650, MG   ASMG iter:  8,  err = 9.1080e-09,  time = 0.22 s
#dof:   7107,  #nnz:  52906, MG   ASMG iter:  8,  err = 3.8331e-09,  time = 0.16 s
#dof:  28547,  #nnz: 216362, MG   ASMG iter:  8,  err = 3.8786e-09,  time = 0.77 s
#dof: 114435,  #nnz: 875050, MG   ASMG iter:  8,  err = 3.6838e-09,  time =  2.9 s
Table: Error
 #Dof        h      |u_I-u_h|_1    ||u-u_h||   ||u_I-u_h||_{max}

  1891   6.25e-02   1.80268e-01   2.14012e-02   3.65459e-02
  7363   3.12e-02   6.32875e-02   5.32944e-03   1.00022e-02
 29059   1.56e-02   2.21618e-02   1.32998e-03   2.61407e-03
115459   7.81e-03   7.77049e-03   3.32204e-04   6.68059e-04

 #Dof        h      ||p_I-p_h||    ||p-p_h||   

  1891   6.25e-02   1.01517e+00   7.11564e-01
  7363   3.12e-02   3.27813e-01   2.20345e-01
 29059   1.56e-02   1.08542e-01   7.02700e-02
115459   7.81e-03   3.68321e-02   2.31884e-02

Table: CPU time
 #Dof    Assemble     Solve      Error      Mesh    

  1891   1.00e-02   2.20e-01   3.00e-02   0.00e+00
  7363   0.00e+00   1.60e-01   2.00e-02   0.00e+00
 29059   6.00e-02   7.70e-01   6.00e-02   0.00e+00
115459   1.40e-01   2.95e+00   1.10e-01   1.00e-02

</pre>
</div>
</div>

<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

<div class="output_area">

    <div class="prompt"></div>




<div class="output_png output_subarea ">
<img src=""
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Conclusion">Conclusion<a class="anchor-link" href="#Conclusion">&#182;</a></h2><p>Optimal second order convergence of velocity in L2 and maximum norm are observed and pressure is observed. For velocity, superconvergence (1.5 order) between nodal interpolate uI and uh is observed. For pressure, the order is also 1.5.</p>
<p>Multigrid solver based on ASMG converges uniformly. Note that due to the bubble function, standard <code>mg</code> solver cannot be used.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div><div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>To-Do: test other boundary conditions.</p>

</div>
</div>
</div>
    </div>
  </div>
</body>

 


</html>
